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We present the phase diagram, in both the microcanonical and the canonical ensemble, of the 
Self-Gravitating-Ring (SGR) model, which describes the motion of equal point masses constrained 
on a ring and subject to 3D gravitational attraction. If the interaction is regularized at short 
distances by the introduction of a softening parameter, a global entropy maximum always exists, 
and thermodynamics is well defined in the mean-field limit. However, ensembles are not equivalent 
and a phase of negative specific heat in the microcanonical ensemble appears in a wide intermediate 
energy region, if the softening parameter is small enough. The phase transition changes from second 
to first order at a tricritical point, whose location is not the same in the two ensembles. All these 
features make of the SGR model the best prototype of a self-gravitating system in one dimension. 
In order to obtain the stable stationary mass distribution, we apply a new iterative method, inspired 
by a previous one used in 2D turbulence, which ensures entropy increase and, hence, convergence 
towards an equilibrium state. 



I. INTRODUCTION 



There are many objects in our universe whose behavior can be understood considering only the gravitational 
interaction. Examples are globular clusters, galaxies, clusters of galaxies, molecular clouds 1]. Different theoretical 
approaches have been proposed to explain the peculiar statistical properties of self-gravitating systems. The main 
difficulty is that these systems cannot approach statistical equilibrium because of the short-distance divergence of the 
potential and of the evaporation at the boundaries. Even if one puts the system in a box with adiabatic walls, thus 
eliminating evaporation, still gravity causes the well-known phenomenon of gravothermal catastrophe |^|^|^. The 
introduction of a small-scale softening of the interaction potential avoids such a catastrophe, so that self-gravitating 
systems can approach the final (thermal) equilibrium state. However, such a state may have a negative specific heat. 
Moreover, a first-order phase transition from the high energy gas phase to the low energy clustered phase appears 

Direct studies of the full three-dimensional A^-body gravitational dynamics are particularly heavy ^ and even 
special purpose machines have been built to this aim Q. Therefore, lower dimensional models have been introduced 
to describe gravitational systems with additional symmetries. For instance, the gravitational sheet model describing 
the motion of infinite planar mass distributions perpendicularly to their surface has been considered (7j. Although 
this model shows interesting behaviors 8, 9], the specific heat is always positive and no phase transition is present. 

Recently, another one-dimensional model has been introduced [Tol | where particle motion is confined on a ring, but 
the interaction is the true Newtonian 3D one. At short distances, the potential is regularized, so that the particles 
do not interact. This model has been called the Self-Gravitating Ring model (SGR) and will be the subject of the 
study discussed in this paper. It has been shown in numerical simulations jlCL] , that this model maintains the peculiar 
features of the 3D Newtonian potential, showing a negative specific heat phase and a phase transition if the softening 
parameter is small enough. Moreover, for large softening, this model reduces to the Hamiltonian Mean-Field model 
(HMF) [ni, which has been recently extensively studied as a prototype system with long-range interactions. This 
latter model, however, although it displays a second order phase transition, does not have a negative specific heat 
phase at equilibrium. 

In this paper, we derive the equilibrium thermodynamics of the SGR model both in the canonical and in the 
microcanonical ensemble. For all non vanishing softening parameter values, this model has a thermal equilibrium 
state. If the softening parameter is small enough, the model shows ensemble inequivalence IT^ with a negative 
specific heat phase in the microcanonical ensemble and a first order phase transition. Therefore, the SGR model 
displays several features of the true 3D Newtonian interaction, and can serve as a better prototype of self-gravitating 
systems in one dimension than all previously introduced models. 

The paper is organized as follows. In Sec. ^ we briefly introduce the SGR model and we discuss the essential 
features of previous numerical simulations 10]. In Sec. IIIII we show the general scheme for deriving all stationary 
density distributions which maximize Boltzniann-Gibbs entropy at fixed total energy and mass. Section llVI presents 
a new iterative method which ensures entropy increase and leads in a unique way towards the stable equilibrium 
single particle distribution function. The method is inspired by a similar one used to compute entropy maxima in 2D 
turbulence 14]. In Sec.0 we describe in full detail how to implement the iterative algorithm in a numerical scheme. 



Figure 1: Self-Gravitating Ring model with a fixed unitary radius. Particles are constrained to move on a ring and therefore 
their location is specified by the angles measured with respect to a fixed direction. Each pair of particles at di and 9j interacts 
through the inverse-square three-dimensional gravitational force. The distance is measured by the chord, as shown in the figure. 



In Sec. IVII we calculate the thermodynamic quantities of the SGR model using the iterative method. We also show 
how, reducing the softening parameter, one enters into a region of ensemble inequivalence, where a tricritical point 
exists which is not the same in the two ensembles [l5j . Finally, in Sec. lVIll we discuss the dynamical evolution of the 
SGR model, emphasizing the properties of relaxation to equilibrium. 



II. SELF-GRAVITATING RING MODEL 

In this section, we briefly present the Self-Gravitating Ring (SGR) model ^^l- this model, particle motion is 
constrained on a ring and particles interact via a true 3D Newtonian potential (Fig. ^ . 
The Hamiltonian of the SGR model is 



H 



1 ^ 1 



2^"' 2N 

1=1 i,j 



V2 1 - cos(0i - Oj) + £ 

where e is the softening parameter, which is introduced, as usual, in order to avoid the divergence of the potential at 
short distances. 

Taking the large e limit, the potential becomes 



2e 



-0{e-^), (3) 



which is the one of the Hamiltonian Mean-Field (HMF) model • It is well known that the HMF model has a 
second order phase transition, separating a low energy phase, where the particles form a single cluster, from a high 
energy gas phase where kinetic energy dominates and the particles are homogeneously distributed on the circle. One 
usually draws the so-called caloric curve, where temperature, given by twice the averaged kinetic energy per particle 
T = fi~^ = 2 {K)/N, is plotted against the total energy per particle U = H/N. In a situation close to that of the 
HMF model, e.g. for e ~ 10, the caloric curve determined from microcanonical numerical simulations is reported 
in Fig. |2Ia). In the homogeneous phase U > Uc{e), the caloric curve is almost linear, while in the clustered phase 
U < Uc{s), it is bent downward. Nonetheless, temperature always grows with energy and one does not observe any 
negative specific heat energy range. However, as it happens for 3D Newtonian gravity simulations j^, when one 
reduces the softening parameter, a negative specific heat phase develops. For instance, in Fig.|21[b), we show two cases 
at small e where three phases can be identified ^3 • 

• a low-energy clustered phase for U < Utopis), where Utop is defined as the energy at which dT/dU = 0. 



• an intermediate-energy phase, Utopis) < U < Uds), with negative specific heat. 
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Figure 2: Caloric curves of the Self Gravitating Ring (SGR) model obtained from numerical simulations of HamiltonianQ Panel 
(a) refers to the softening parameter value e = 10, for which a second order phase transition appears at Uc- No backbending of 
the caloric curve, indicating a negative specific heat, is present. Simulations were performed for A'^ = 100. Panel (b) shows the 
caloric curves for two different values of the softening parameter, £i = 1.010"^ and £2 = 2.510~^, and = 100. The transition 
is here first order in the microcanonical ensemble (see Sec. I VII for a discussion). The two transition energies Uc{si) and Uc{e2) 
are pretty close, suggesting a slow variation of the critical energy with the softening parameter e. On the contrary, Utop{si) is 
significantly smaller than Utop{£2), indicating that this characteristic energy value diminishes with e. A negative specific heat 
phase appears for Utop <U < Uc, and expands as the softening parameter is reduced. 



• a high-energy gaseous phase for Uc{e) < U. 

The clustered phase is created by the presence of softening e, without which the particles would fall into the zero 
distance singularity. In the gas phase, the particles are hardly affected by the potential and behave as almost free 
particles. The intermediate phase is expected to show the characters of gravity, persisting and even widening in the 
e ^ limit. 

In the following, several of these features will be given a theoretical explanation and we will detail the analysis of 
the nature of the phase transition (first or second order) when e is varied. 
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III. STATIONARY DENSITY DISTRIBUTION 



In the mean- field limit (N oo with fixed length [l^'). one can introduce the single particle distribution function 
f{p, 9) such that /(p, 0)dp(i9 is the fraction of particles in the domain [9, 9 + d9] [p, p + dp] . In terms of /, the potential 
energy can be written as 

Ep[f] = ^Jdedct>dp dp' f{p, 9)V,{9 ~ c^)f{p\ 4>) (4) 
= i / d(lyd9 p{9)p{cj,)V,{e-cj,) (5) 



EK[f] = l I dedpp'f{p,9) (7) 



where 

Pi0)^ I dpfip,9) (6) 

is the mass density. The kinetic energy is 

Ez^lf] = ^ 
and the total energy 

E[f]^EK[f]+Ep[f]. (8) 
The equilibrium distribution in the microcanonical ensemble is determined by maximizing entropy 

S[f] = - Jd9dpf\ogf (9) 

under the constraints of fixed total energy, momentum and mass. In the following, we fix the total energy E[f] = U, 
the total mass 



M[f] = J pd9 = l (10) 

and the total momentum 

P[f] = J Pfip,9)d9dp = 0. (11) 
A necessary condition to get an entropy maximum is to require that the free energy 

F[f] = S[f]- m.f] ^aj fdpde- jp[f] , (12) 
where a, (3 and 7 are Lagrange multipliers, is stationary 

'^^[^1 = _ log/ - 1 - /3 f ^ + W{9)] ^a-jp = 0, (13) 



Sf "V2 
where W{9) is defined as 

w{0)= r\{(t>)v,{0 ~ ^)d(t>. (14) 

J — 7T 

Since p[f] — 0, the Lagrange multiplier 7 vanishes. From Eq. H13() . the normalized stationary distribution function 
can be written as 



f{p,9)=Ae^p 



y + m^) 



(15) 
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where A = exp(— 1 — a) is the normaUzation constant and the mass density is given by 

p(0)=Ie-^^(»), (16) 
where A = Ay/2TT/p. When Eq. (|14|l and H16|) are combined, we obtain the consistency equations 

W{e) = A r e~'^^^^^V,{9 - <j>)d(l), (17) 



and the equihbrium density equation 



p{6) = Aexp 



PA p{cb)V,{e ^ cj,)dcb 



(18) 



which are solved numerically in the following. Once the stationary mass distributions p and the function W are 
obtained for each value of e, the full single particle distribution function f{9,p) is derived from Eq. H15() . The 
potential energy and the kinetic energy are determined by Eq. (jSJ and Eq. (O respectively, allowing to draw the 
caloric curve by plotting T = f3~^ — 2Ek against the total energy U — Ek + Ep. 



IV. AN ITERATIVE METHOD TO SOLVE THE EQUILIBRIUM DENSITY EQUATION 



The inverse temperature /3 can be expressed in terms of the energy 

''+7r /'+7r 



TT J — TT 



(19) 



Once an initial density distribution p(^,(&) is chosen, one can determine an initial inverse temperature /3o using Eq. H19|l . 
and then solve the consistency equation H18|l iteratively (as done for instance in Ref. 0]). However, we will follow here 
a different iterative method, which ensures entropy increase and, hence, convergence of the algorithm. The method 
is inspired by a similar one used by Turkington and Whittaker [ij to compute entropy maxima for two dimensional 
turbulence. 

The functional to maximize is strictly concave and we must fix both a linear constraint M[/] — 1 and a 
nonlinear one E\j\ = U . It is this latter nonlinear constraint which makes the variational problem more difficult than 
usual. The trick to solve this nonlinear problem consists in considering a linearization of the energy constraint around 
the distribution function resulting from the previous step in the iterative process. 

One begins with the normalized distribution obtained at the fcth step of the algorithm. From that, one can 
compute the mass density pk and the average potential VFfc. 



pm = J dpfkip,9) 
Wk{9) = r dcf,pk{^)V,{9-cl>). 

J —TT 

The distribution at the next step fk-\-i wiU be then determined by solving the fohowing variational problem 

{f-fk)dpde<u], 
J "J h J 

where the functional derivative of the energy is 

Wk {9) . 



max|^[/] I M[f] = l,E[fk]+ J ^ 



SE 



2 



(20) 
(21) 

(22) 
(23) 



This variational problem has a unique solution fk+i, since it corresponds to the maximization of a strictly concave 
functional with linear constraints. 

This iterative process ensures convergence of the entropy. Let us prove it. By using a generalization of the Lagrange 
multiplier rule for our inequality constrained variational problem [iM Il9j 



5S_ 
'5] 



fk + l 



SE 

= ttfe+i + fJk+i -J7 



(24) 
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with the additional requirement 



+1 



SE 



ifk+i- fk)dpd9~U 



= 0, 



(25) 



where Pk+i > is the multipher associated with the energy constraint and a^+i, the one associated with mass 
conservation. When solving Eq. H25(l . we have either /3fc+i — 0, which removes the energy constraint, or Pk+i > 0, 
and an equality for the linearized energy constraint. 

In order to prove convergence of the entropy, let us first prove that the energy functional E [f] is concave. Since 
the kinetic part is linear in /, the second variation of E[f] is 



6^E 



5p(e)5p{(j>)v,(e-<j>) 



E ^-.'^ \^p'^ 



(26) 
(27) 



where the second equality is obtained using the Fourier series expansion for both the mass density variation 5p and 
the potential Ve 



5pk = 
V.k = 



1 

2^ 
1 

2^ 



dip exp [ikip) dp (ip) 
dip exp (ikip) Ve (ip) . 



(28) 
(29) 



Since is even in the argument (j), Ve,k is a real number. Moreover, since < and Ve (<p) is strictly increasing for 
< (/J < TT, it is easy to prove that for any fc, V^^k is strictly negative. Hence, from formula (|27|l the second variation 
of the energy functional is negative and this functional is strictly concave. 

On the other hand, the entropy is strictly concave. We have 



S[f + 5J]<S[J] 



d9dpf- 



Sf-l I dOdp 



f 



f 



(30) 



where in the derivation we have used ln(l + a;) > a; — x"^ /2 for x > —1. Applying this latter inequality with / = fk+i 
and Sf~fk — fk+i, and using both condition l|24|l and 1)25(1 . we obtain 



S [fk+i 



S[fk]>Pk+i{U-E[M) + - I d9dp 



(/fe+i - /fe)^ 



fk 



(31) 



where the term proportional to a^+i vanishes because of mass conservation. 
We will now use the concavity of the energy functional E [/]. For fc > 1, 



E[fk]<E[fk-i]+ J 



5E_ 
Sf 



ifk- fk~i)dpd0. 



As Pk+i > and E [f^] < U, directly from the variational problem 



S[fk+i]-S[h]> 



(/fc+i - fkY 

fk 



Eq. (j2U implies that 



dedp > 0. 



(32) 



(33) 



Hence, the entropy has to increase for all iterates after the second. Since the entropy is bounded from above, it has to 
converge. Using Eqs. H31() and H33() . one derives that the energy E [fk] converges to U from below. Moreover, assuming 
that fk converges toward /, one can prove the convergence of the multipliers to limit values a and /3 > 0, which implies 
that / verifies Eq. ((24|l for equilibrium states. Although mathematically one cannot prove the convergence of /, in 
all practical cases we will analyze, it appears to be verified. For a more thorough discussion of the convergence in the 
similar case of the Euler equation, see Sec. IV in Rcf. [T^ . 
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V. IMPLEMENTATION OF THE ALGORITHM 



We describe in this section the practical implementation of an algorithm which allows the calculation of the stable 
distribution, using the method described in the previous section. 
From H24|) . we obtain 

(34) 



(35) 



fk+i = Ak+i exp y-(3k+i + 
where Ak+i = exp(— a^+i — 1) and f3k+i are unknown at this stage. Using H20() . we get 



Pk+i{0) = Ak+i e 



where Ak+i = Ak+iy/^n/ (3k^i. This equation allows to compute Wk+i{0) from Eq. and 



k+l 



E[fi 



fc+ij 



Pk+i{9)Wk+i{0)d0. 



(36) 



Then the multipliers ak+i and Pk+i must be computed from Eqs (|10(l and (|25() and, from these, one gets Ak+i- In 
order to compute numerically these Lagrange parameters, let us define the Lagrangian 



Lfc [/](/?, a) = -5[/]+/3 
From this, one further defines 



if-fk)dpd9-U 



a(M[f] - 1). 



L^(/3,a) =inf{Lfc[/](/3,a)}. 



(37) 



(38) 



One can prove on a general ground '19^ that is concave and that ak+i and Pk+i are the unique maxima of L^. 
Using condition (|34|l for the extrema of Lfc[/](/3, a), we can compute L^. We obtain, using for practical reasons the 
variable A instead of a, 



1 



Lli(3,A) = logA+-logP-f3[U + Ek- — 



1 



A 



Necessary conditions for the concave function to be maximal are 

'd(?e-'3^'=W =0, 



dI4 
dA 



1 

I 
1 



— TT 

1 



U-Ek + A 



dp 2/3 2(3k 
Substituting Eq. H4U|) into Eq. H41|l . one gets the condition 



(39) 

(40) 
(41) 



1 



2/3fe+i 213k 



T^ + U + Ek 



dO H/fe(6')e-^'=+i^'=W 



0, 



(42) 



which, since is concave, has a unique solution. Numerically, the solution pk+i is found by using a Newton algorithm 
for Eq. (|42|) . Then, from Eq. H4U|) . we get Ak+i- Finally, we can calculate the new density distribution from Eq. 135|l . 



VI. DISCUSSION OF THE RESULTS 



Using the iterative method described in the previous section, we are able to derive the stable mass density p{0) 
solution of Eq. H18|l and, from that, all thermodynamic functions in the microcanonical ensemble. In the first part of 
this section, we will show the numerical solution obtained for p{9), and its dependence on energy for a small value 
of the softening parameter e. In the second part, we will discuss the phase diagram of the SGR model, both in the 
microcanonical and in the canonical ensemble, when e is varied. 
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A. Mass density, entropy and caloric curves 

For energies above a certain critical value Uc{s), the stable mass density solution is uniform. In this case, one can 
compute the entropy from Eq. ^ 



5 = i(31og(27r) + l-log/3), 



(43) 



and the inverse temperature from Eq. 



f3 = {2U~2Epy 



(44) 



where 



Ep — 



1 1 



2 (27r)2 
1 



7rV2 \/2" 



IC 



d9d(l)V, 
2 



(45) 
(46) 



Ep, because this latter 



where K, is the complete elliptic integral of the first kind IC{x) = /J^^^ d9/ \/l — x sin^ 6. 

Remark that Eq. H44|l implies that the homogenous state cannot be continued below Uhom 
energy corresponds to zero temperature. 

For U < Ucie), the stable mass distribution must be determined numerically. We have checked in this case, that a 
direct iterative method of solution of the consistency Eqs (|17() and H18() does not always converge. On the contrary, 
the novel algorithm presented in Sec. ensures convergence as shown in Fig. 01 for the entropy. 




10 20 30 40 50 60 70 80 90 100 



Figure 3; Convergence of the entropy using the algorithm of Sec.EIfor £ = 10 ^ ,U 



In Fig. 21 we show both entropy and temperature T = as a function of energy U. The most striking feature 
is the presence of a negative specific heat region for Utop < U < Uc- For Uiow < U < Uhigh, the entropy does not 
coincide with its convex envelope. Hence, microcanonical and canonical ensemble do not give the same predictions. 
Indeed, the main peculiarity of the microcanonical ensemble is that macroscopic states within this interval are stable, 
while they would be either metastable or unstable in the canonical ensemble. The mass density is uniform above 
Uc, while, below this value, it is localized. The appropriate order parameter to characterize this localization is the 
"magnetization" 

B = r dee'%{e), (47) 

which vanishes if the mass distribution is uniform while it reaches the value B — I when the mass is concentrated in 
only one point. Intermediate degrees of localization give intermediate values of B. The "magnetization" is plotted in 
Fig.|31as a function of U. It is a decreasing function of U, up to Uc, where it has a jump to the limiting value 0. Hence, 
we have here a first order microcanonical phase transition. The first order nature of the phase transition is confirmed 
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Figure 4: Temperature (panel (a)) and entropy (panel (b)) versus energy U for the softening parameter value e = 10~^. Four 
values of the energy, indicated by the short-dashed vertical lines, can be identified from this picture: Uiom — —93 and Uhigh — 6 
bound from below and above the region of inequivalence of ensembles. f7c — is the transition energy in the microcanonical 
ensemble. Utop — —66 limits from below the negative specific heat region, where temperature decreases as energy increases. 
Tcan — 15, represented with a dashed line in panel (a), is the canonical transition temperature and corresponds to the inverse 
slope of the entropy, both at Uio-w and Uhigh, as represented by the straight dashed line in panel (b). The full lines represent 
the analytical solutions of the temperature and of entropy in the uniform case (see formulas (1431 and 1441 ) . They are extended 
slightly below Uc, in the metastable phase, in order to identify them. The insets in panels (a) and (b) show a zoom of the 
temperature and of the entropy around Uc, revealing a temperature jump at Uc and different slopes of the entropy above and 
below Uc, which emphasizes the first order nature of the phase transition. 



zooming the entropy around Uc (see the inset in panel (b) of Fig. 2)). This reveals that this first order phase transition 
is of the convex-concave type (see Ref. ) . The canonical ensemble is obtained by taking the convex envelope of the 
microcanonical entropy. The transition is first order in the canonical ensemble and the transition temperature Tcan is 
given by the inverse slope of the entropy at Uiow and Uhigh ■ No canonical macrostate is present in the energy range 

[Ulow, Uhigh]- 

A typical localized mass density distribution is shown in Fig.El It corresponds to an energy where the specific heat 
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Figure 5: "Magnetization" B versus energy U for e = 10 ^, which emphasizes the microcanonical first order phase transition 
at ?7c — by showing a jump in the order parameter. 
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Figure 6: A typical mass density distribution p{d) for e = 10~^ and U — —20.0, in the negative specific heat region. 

The first order phase transition is associated with the existence of metastable states. Using a continuation method, 
we have been able to compute them. Their entropy is represented in Fig. \7\ around the transition energy Uc for the 
specific case e ~ 10^^. The inhomogeneous metastable state turns out to exist for Uc < U < Um with Um — 0.16, 
while the homogeneous metastable state exists for Ufiom I^U < Uc, with Uhom = Ep{e ^ 10^^) ~ —1.19. 

B. Behaviour as the softening parameter e is varied 

Let us first examine a situation where the softening parameter is much larger than previously, e = 10^^. In 
the microcanonical ensemble. Fig. |S1 shows that a concavity change still occurs at Utojy — —0.8, and that a phase 
transition exists sd U — Uc — —0.3. However, the temperature being now a continuous function of the energy but 
with discontinuous derivative at Uc, the phase transition is of second order, and is associated with the symmetry 
breaking of the order parameter. The caloric curve shows that this second order phase transition is of the convex- 
concave type. As it is necessary for this type of microcanonical second order phase transition jl3(|, we observe a 
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u 

Figure 7: Both the high energy branch of the entropy versus energy curve, corresponding to the homogeneous solution (solid 
line) and the low energy branch of the inhomogeneous solution (dashed line) are represented in this plot for e — 10~^. The two 
branches cross at Uc — 0. The continuation of the homogeneous branch into the low energy region is bounded from below by 
Uhom — —1.19, indicated by a vertical dashed line. The inhomogeneous branch continues into the high energy phase and ends 
at Uin — 0.16, again indicated by a vertical dashed line. 



positive specific heat jump at the transition point. 

What we find suggests that between e = 10~^ and e = 10^^, there is an intermediate value of e where a micro- 
canonical tricritical point is present. This point is signalled by two properties: 

• The caloric curve assumes a negative infinite slope as U tends to Uc from below. 

• The upper energy of the metastable inhomogeneous phase Uin collapses onto Uc from above, while still a 
continuation of the homogeneous phase below Uc exists as an unstable phase. 

In Fig. 1101 we have represented the e-dependence of the critical energy Uc and of the energy bounds Uin and Uhom- 
At the microcanonical tricritical point, ~ 10~^, the end point for the existence of the inhomogeneous metastable 
phase joins the critical line. This is a generic feature of tricritical points with symmetry breaking (see Fig. 6 of 
Ref. ;i3]). 

To locate the tricritical point in the canonical ensemble, one has to look for the e- value at which the two curves Uiow 
and Uhigh merge (see Fig. lllll . An approximate estimate of this value is e9p ~ 10~^. At the canonical tricritical point, 
also Utop merges with the above curves, indicating the disappearance of the negative specific heat region. We thus 
note that ensemble inequivalence disappears at the canonical tricritical point by the disappearance of the inflexion 
point at Utop in the entropy curve. As it may be checked in Table I of Ref. '13^, this is the only way in which ensemble 
inequivalence can disappear when associated with a tricritical point. 

Summarizing, the important changes of the phase diagram of the SGR model when e is varied are due to the 
existence of microcanonical and canonical tricritical points. For e < e^, there is an energy range with ensemble 
inequivalence. These features were already observed in Refs. [l5ll20| . 



VII. RELAXATION TO EQUILIBRIUM 



We have first checked numerically if the equilibrium density profile is ever attained in direct iV-body simulations 
of Hamiltonian In Fig.^l we compare the result of a numerical simulation with the equilibrium density profile 
obtained by the iterative method. The agreement is good in the center of the mass distribution, while the tails are 
still affected by strong flnite N fluctuations. In this case, the total energy U = —20 is in the region of negative 
speciflc heat and is well conserved using a sixth order symplectic integrator plj . Simulations were performed using 
GRAPE-5, a special purpose computer for gravitational force 0]. 

However, it is well known that systems with long-range interactions display a very slow relaxation to equilibrium |22j| . 
Hence, we expect that similar features will be also exhibited by the SGR model. For instance, we can consider a 
"cold start", where the particles are initially homogeneously distributed on an arch {9 £ [OmimOmax]) with zero 
kinetic energy. Usually, in gravitational simulations, one looks at the evolution of the virial ratio \2K/V\, which is 
here initially zero. The plot of the time-dependence of this parameter is shown in Fig. 1131 for the same parameter 
values used previously. One clearly observes that the system relaxes to a "quasi-equilibrium" state, where the virial 
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Figure 8: Panel (a): The caloric curve (triangles) for e = 10~^. The dashed vertical lines indicate, from left to right, 
Uiow — —1.98, Utop — —1.3, Uc — —0.32 and Uuigh — —0.225. The homogeneous phase curve, known analytically, is shown by 
the continuous line. The main difference with respect to Fig. 0] is that now there is not a temperature jump at Uc- The phase 
transition is second order in the microcanonical ensemble, while it is still first order in the canonical ensemble, at Tcan — 0.8. 
Panel (b): Entropy versus energy (triangles) for e = 10~^. The entropy curve corresponding to the inhomogeneous distribution 
smoothly connects with the one of the homogeneous distribution (solid line) . The oblique straight dashed line is tangent to the 
entropy at Uiow and Uhigh, which delimit the energy region of ensemble inequivalence. 



ratio fluctuates around a value which differs from the equihbrium one, computed analytically. While previously, for 
the mass positions (see Fig. I12II . the relaxation was observed on a short time scale, we show here that a quantity 
related also to velocities does not display a relaxation on the same timescale. From previous experiences with similar 
cases li^l, one expects that the relaxation should occur on a time scale of the order of a power of N. 

Even slower is the relaxation when local maxima of the entropy exist. This happens around the critical energy 
Uc in the case of a first order phase transition, e.g. for e = 10~^. Figs 1141 show the relaxation to different values 
of two relevant quantities, the temperature and the "magnetization" , when the system is initialized either with the 
particles concentrated on a small arch, or on a larger one. When the system is "close" to the local entropy maximum 
corresponding to the clustered state, it converges to it pretty fast. The contrary happens when the particles are more 
homogeneously distributed, and then the system converges to the homogenous state. Indeed, between the two states 
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Figure 9: "Magnetization" B versus energy U for e = 10 which emphasizes the microcanonical second order phase transition 
at Uc, because the order parameter vanishes continuously. 
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Figure 10: The dash-dotted line represents Uinie), the dashed line Uhom ~ Ep{e), the filled circles the first order microcanonical 
phase transition energy and the open circles the second order one. At the microcanonical tricritical point ~ 10~^, the phase 
transition changes from first order to second order and, at the same time, the inhomogeneous metastable solution disappears. 



there is an entropy barrier which has been found to grow as exp (iV) for systems with long-range interactions [23t |2 



VIII. CONCLUSIONS AND PERSPECTIVES 



We have fully characterized from the thermodynamic point of view a one-dimensional model of self-gravitating 
particles moving on a ring [loj . which is the simplest prototype of the full 3D self-gravitating system. Solving the 
equilibrium density equation by a new iterative method, whose convergence is assured by entropy increase, allows 
to derive the full phase diagram of the model both in the microcanonical and the canonical ensemble. When the 
softening parameter is sufficiently small, a negative specific heat region appears in the microcanonical ensemble, in 
coincidence with the phase transition becoming first order in the canonical ensemble. Further lowering the softening 
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Figure 11: e-dependence of Uiom (dashed line), Uhigh (solid line) and Utop (dash-dotted line). The canonical tricritical point is 
located at ~ 10~^ where the three curves merge. At this softening parameter value also the negative specific heat disappears 
in the microcanonical ensemble, while the transition becomes second order in the canonical ensemble. In the figure, we also 
show, with a long-dashed line, the theoretical estimate of Ut^p ~ — l/(4v'2e) obtained in Ref. |lO^ . 
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Figure 12: Comparison of the mass density profile obtained by the iterative method (solid line) with the result of numerical 
simulations (plus signs) with A'' = 4000 and e = 10~^. Parameter values are the same of Fig|S] The inset is a zoom of the 
center of the profile. 




parameter, the transition becomes first order in the microcanonical ensemble and a temperature jump appears at the 
transition energy. The microcanonical and canonical tricritical points do not coincide [l5j. 

Dynamically, we have performed numerical experiments which show that relaxation to equilibrium can be extremely 
slow. They reveal also the presence of quasi-equilibrium states, which are ubiquitous in systems with long-range 
interactions [l^ . These states could be further characterized considering a Vlasov equation approach as done for the 
HMF model in Ref. [22|. Moreover, in the first order microcanonical transition region a strong metastability appears 
and, at a given energy, the system can relax towards different thermodynamic states. 

Preliminary studies of velocity probability distributions in this model have been performed in Ref. 10]. A similar 
analysis has been more recently done for the full 3D self-gravitating system 12^. In both models, non gaussians tails 
show up in several energy regions. Among future directions of study of the SGR model, we think that deriving a 
theoretical framework to undertand these large tails would be of particular interest. To this aim, especially useful 
could be the methods developed to understand single particle diffusion in the HMF model [2^ |27| . 
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Figure 13: Time evolution of the virial ratio \2K/V\ of the SGR model for e = 10~^ and N = 4000. Initially, the particles are 
homogeneously distributed in the interval [0, 27r/75] with zero kinetic energy. The virial ratio oscillates asymptotically around 
the value 0.49, which differs significantly from the equilibrium value 0.55 indicated by the dashed horizontal line. The initial 
virial ratio is zero, although this time region is not visible in the figure. 
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Figure 14: Relaxation to difi'erent maximum entropy states in the SGR model for U ~ 0, e — 10^"' and N — 10'^. Panel (a) 
shows the relaxation of the temperature either to the inhomogeneous state value (horizontal dotted line), or to the homogeneous 
one (horizontal dash-dotted line), depending whether the particles are initially distributed on a smaller arch 9 £ [0, 7r/50] (solid 
line) or a larger arch 9 £ [0, tt/S] (dashed line). In both cases, the velocity distribution is initially a 'Svater bag". Panel (b) 
shows the same for the "magnetization" . 
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